Critical temperature of interacting Bose gases in two and three dimensions. 
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We calculate the superfluid transition temperature of homogeneous interacting Bose gases in 
three and two spatial dimensions using large-scale Path Integral Monte Carlo simulations (with up 
to N = 10 5 particles). In 3D we investigate the limits of the universal critical behavior in terms 
of the scattering length alone by using different models for the interatomic potential. We find that 
this type of universality sets in at small values of the gas parameter na 3 < 10 -4 . This value is 
different from the estimate na 3 < 10 -6 for the validity of the asymptotic expansion in the limit of 
vanishing na 3 . In 2D we study the Berezinskii-Kosterlitz-Thouless transition of a gas with hard-core 
interactions. For this system we find good agreement with the classical lattice \4>\ 4 model up to very 
large densities. We also explain the origin of the existing discrepancy between previous studies of 
the same problem. 

PACS numbers: 



The theoretical determination of the superfluid transi- 
tion temperature in homogeneous, interacting Bose sys- 
tems is a fine example of a many-body problem that can 
be quantitatively addressed only by "exact" numerical 
techniques. This fact is well understood in the case of 
strongly interacting quantum fluids, such as liquid 4 He 
where the critical temperature in bulk [l[ as well as in 
two-dimensional configurations was calculated us- 

ing path integral Monte Carlo (PIMC) simulations, but 
at first glance is surprising in the case of dilute gases. 
However, in three dimensions (3D) the presence of any 
finite interaction changes the universality class of the 
transition from the Gaussian complex-field model, cor- 
responding to the ideal gas Bose-Einstein condensation 
(BEC) temperature T c °, to that of the XY model. Thus, 
the critical temperature T c can not be obtained from Xj? 
perturbatively [4|. In two dimensions (2D) the superfluid 
transition, which belongs to the Berezinskii-Kosterlitz- 
Thouless (BKT) universality class [1, Q , is induced by 
interaction effects and there is no unperturbed critical 
temperature to start with. 

In a 3D weakly repulsive gas the critical temperature 
shift is fixed by the s-wave scattering length a (a > 0), 
which characterizes interatomic interactions at low tem- 
perature 0, 0, s s m in m m , 



T, = T° 



l + c(an^ 3 ) 



(1) 



Here n is the gas density and T® = 
(27rfi 2 /mfc B )[n/C(3/2)] 2 / 3 with m the particle mass 
and C(3/2) — 2.612. The numerical coefficient c in 
Eq. (jTJ) was calculated in Refs. 0, EH by solving the 
effective 3D classical \ip\ 4 model using lattice Monte 
Carlo simulations. The reported value is c = 1.29 ±0.05. 
The same classical model was employed in Ref. [ll| to 



calculate the higher-order logarithmic correction to (JTJ) . 

Continuous-space studies of a gas of hard spheres, 
based on the conventional PIMC algorithm [l4|, were 
carried out in Refs. 0, EH. Both calculations suffered 
from two shortcomings: first, the number of particles in 
the simulation was only few hundreds making the extrap- 
olation to the thermodynamic limit difficult; second, the 
algorithm is known to be inefficient for simulations of the 
superfluid density. The results of Ref. [T3| for T c agree 
with the asymptotic law ([T]) at small densities, in con- 
trast to the significantly lower values obtained in Ref. 0] 
(see Fig. [T]). Moreover, there is strong (about ten stan- 
dard deviations!) discrepancy between Refs. 0] and [T3| 
at higher densities calling for further investigation of the 
problem. 

In 2D the BKT transition temperature of a weakly 
interacting gas is written in the form (TBI 



Tbkt = 



1 



mkr 



ln(£/4^)+lnln(l/n 2I)a 2 D ) » ( 2 ) 

where ri2Da2D ^ s * ne corresponding gas parameter, and 
a2D is the 2D scattering length. The numerical coefficient 
£ was calculated in Ref. [I6( from lattice Monte Carlo 
simulations of the classical \ip\ 4 model, similarly to the 
3D case, yielding the value £ = 380 ± 3. 

An important question concerns the range of applica- 
bility of Eqs. H])-© since they were derived for the limit 
of vanishingly small interaction strength. More broadly, 
one is interested in knowing up to what value of the gas 
parameter it is possible to express T c as a function of 
na 3 alone and ignore the dependence on the interaction 
potential details. These questions are particularly rele- 
vant for the 2D case where experimental determinations 
of the critical tem per ature in trapped configurations have 
become available 17, 1811. 
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In this Letter we report PIMC results for the superfluid 
transition temperature in interacting Bose gases both in 
3D and in 2D. We carry out large-scale simulations of 
homogeneous systems with up to N = 10 5 particles. The 
simulations are based on the worm algorithm, recently 
extended to continous-space systems, which allows for a 
reliable and efficient calculation of the superfluid den- 
sity [![. With new data the extrapolation of the critical 
temperature to the thermodynamic limit can be done 
with the level of accuracy that was unreachable in pre- 
vious attempts. In 3D we determine the dependence of 
T c on the gas parameter na 3 , reporting good agreement 
with the expansion ([1]) in the dilute regime and signif- 
icant deviations from previous studies [3, EH at higher 
densities. Furthermore, we carry out simulations with 
both a hard- and a soft-sphere interatomic potential to 
investigate the universal behavior in terms of the scat- 
tering length. In 2D we calculate Tbkt for a hard-disk 
ffelS 0jS 0j function of the interaction strength, finding re- 
sults in excellent agreement with the prediction (|2|) up to 
regimes of surprisingly high density. 
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FIG. 1: (color online). Critical temperature in 3D as a func- 
tion of the gas parameter na 3 . The symbols labeled by PRA04 
correspond to the results of Ref. the ones labeled by 

PRL97 correspond to Ref. 0. The dashed line (green) is the 
expansion JTJ of Ref. [I| and the dotted line (black) is the 
expansion of Ref. [llj including logarithmic corrections. The 
solid line (red) is a guide to the eye. 



We consider a 3D system of N particles with peri- 
odic boundary conditions described by the Hamiltonian 



H = -(^ 2 /2m)Eti V? + Ei<j V(\vi - r,|), where r, 
denotes the coordinates of the i-th particle. Two-body 
interactions are modeled by the following potentials: the 



hard-sphere (HS) potential, V (r 



if 



r < a 

and zero otherwise, and the soft-sphere (SS) potential, 
V ss (r) — Vo (Vq > 0) if r < R and zero otherwise. In 
the HS case the s-wave scattering length a coincides with 
the range of the potential, while in the SS case it is given 
by a = R [l - ta,nh(K Q R Q )/K Ro] with K = ^/V^m/h. 
We always use the range i?o = 5a and adjust Vq to ob- 
tain the desired value of a. We notice that the HS and 
SS model represent two extreme cases of repulsive po- 
tentials: in the HS case the energy is entirely kinetic, 
while for the SS case the Born approximation result 
OB = (m/ti 2 ) J Q V(r)r 2 dr accounts for more than 80% 
of the value of the scattering length. 

In a PIMC simulation one obtains averages of physical 
quantities over a set of stochastically generated configu- 
rations R = (ri, ...,rjy) sampled from a probability dis- 
tribution proportional to the density matrix p(R, R, (3) = 
(R|e~^ ff |R), where (3 — 1/fcsT is the inverse temper- 
ature. The superfluid density p s is obtained from the 
winding number estimator which accounts for long 
permutation cycles of identical particles occurring in the 
system. The calculation of /o(R, R, 0) is based on the 
pair-product decomposition, where the two-body density 
matrix associated with the relative motion of the pair is 
determined exactly both for the HS and the SS poten- 
tial 0,11. The critical temperature T c is determined 
from calculations of the superfluid fraction p s / p (p = mn 
is the total mass density) for systems with increasing par- 
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FIG. 2: (color online). Results for a 3D non-interacting gas. 
Scaled superfluid fraction as a function of temperature for 
different system sizes. The lines are linear fits from Eq. ©. 



ticle number N using the scaling ansatz[21| 

N^p^t, N)/p = f(tN^ v ) = /(0) + f'^tN 1 ^ + ... . 

(3) 

Here, t = (T — T c )/T c is the reduced temperature, v is 
the critical exponent of the correlation length £(t) ~ i - ", 
and f(x) is a universal analytic function which allows for 
a linear expansion around x = 0. 

We first consider the non- interacting case (see Fig. [2]). 
The scaling curves all intersect at the same value of the 
reduced temperature according to Eq. ([3]) with T c /T° = 
1.0005(4) and v = 0.96(3). The value of v is consistent 
with the prediction v — I of the complex Gaussian model. 
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1.058 1.06 1.062 1.064 1.066 1.068 1.07 1.072 1.074 
T/T° 



TABLE I: Transition temperature T c of a 3D hard-sphere gas 
for different values of the gas parameter. The results with the 
label SS correspond to a soft-sphere gas. 
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x 10" 
5 x 1(T 7 
1 x 1(T 6 

i x irr 5 

1 x 1(T 4 
(SS) 1 x 10" 4 



T c /T° 



1.0069(10) 
1.0091(7) 
1.0127(7) 
1.0214(7) 
1.0351(7) 

1.0359(12) 



2 x 10" 3 
(SS) 2 x 10~ 3 

5 x 10~ 3 

1 x 10 -2 
(SS) 1 x HT 2 

5 x 10" 2 



T c /T° 



1.0624(4) 
1.0277(4) 
1.0652(4) 
1.0627(4) 
0.9880(5) 
1.0060(5) 



FIG. 3: (color online). Results for a 3D hard-sphere gas with 
na 3 = 5 x 10~ 3 . Scaled superfluid fraction as a function of 
temperature for different system sizes. The lines are linear 
fits from Eq. @. The inset shows the N dependence of the 
intersection point between lines corresponding to pairs of con- 
secutive system sizes. 



The corresponding results for the hard-sphere gas at the 
density na 3 — 5 x 10 -3 are reported in Fig. [3] In this 
case the intersection point between curves for consecutive 
system sizes clearly drifts towards lower temperatures for 
larger N (see inset of Fig. (U). This effect arises from cor- 
rections to the scaling law [the right hand side of Eq. ([3|) 
has to be multiplied by (1 + CN~ U / 3 + . . . ) with lo m 0.8 
and C of order unity]. A reliable extrapolation to the 
thermodynamic limit using Eq. ([3]) requires large-scale 
systems on order of N > 10 4 . The temperature T c and 
the exponent v are determined from Eq. ^ by consider- 
ing series of data corresponding to systems with A~ suf- 
ficiently large. We find T c /T° = 1.0652(4) (see Table |J) 
and v — 0.70(4), in agreement with the value v = 0.67 
corresponding to the universality class of the XY model 
in three dimensions [22| ■ By using the same procedure 
we calculate T c as a function of the gas parameter na 3 
both for the HS and the SS potential. The results are 
shown in Fig. [1] and are reported in Table fl] 

Except for the largest densities, our results for the HS 
gas are systematically higher than the ones of Ref . and 
lower compared to the ones of Ref. [l]| [23| • With much 
better accuracy for p s and larger system sizes we are in 
the position to explain the origin of discrepancy. It is 
two-fold. First, the number of imaginary time slices used 
in Ref. [l3| was 15, a factor of 3 larger than in Ref. Q • We 
find, by doing simulations with up to 96 slices, that about 
25 slices have to be used to ensure that the correspond- 
ing systematic errors are negligible. Second, Refs. [3, EH 
underestimated error bars by a factor of ten because mul- 
tiple intersections between scaling curves (i) render the 
procedure of locating the intersection point ambiguous, 
and (ii) prevent one from detecting corrections to scaling 
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4: (color online). Results for a 2D hard-disk gas with 
na~ = 0.01. Superfluid density as a function of tempera- 
ture for different system sizes. The solid line (red) shows 
the extrapolation to the thermodynamic limit. The dashed 
line (black) corresponds to the BKT universal jump. In the 
inset we show the dependence of the superfluid density on 
the system size for different temperatures. The dotted lines 
are fits using the Kosterlitz-Thouless renormalization group 
equations. 



and the flow of intersection points with the system size 
(which is substantial between A" = 256 and A~ ~ 20000). 
The discrepancy with Ref. [T3| reduces at very small den- 
sities na 3 < 10 -6 , where finite-size corrections to T c ap- 
pear to be less relevant and there is agreement with the 
expansion ([1]) . The critical temperature T c first increases 
with na 3 , then goes through a maximum and for larger 
values of the gas parameter decreases below the T° value. 
For example, T c /T® = 0.70 in liquid 4 He corresponding 
to the effective gas parameter na 3 ~ 0.21 [2J|. The HS 
gas is finally expected to become a solid at the freezing 
density na 3 ~ 0.25 (24[. Concerning the comparison be- 
tween different model potentials, at na 3 = 10~ 4 we find 
very good agreement between the HS and SS gas, while 
for higher densities deviations start to become evident 
and the SS results are significantly smaller. 
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FIG. 5: (color online). Critical temperature of a 2D gas of 
hard disks as a function of the gas parameter na 2 . The solid 
line is the result ([2j|. 



The simulations of the 2D Bose gas are carried out 
in a similar way, using a hard-core potential with range 
«2D (hard-disk potential) to model the interactions. The 
results for the superfiuid fraction as a function of temper- 
ature are reported in Fig. [5] for the value U2do^ d = 0.01 
of the gas parameter. The extrapolation to the thermo- 
dynamic limit is carried out by fitting numerical data at 
finite N to the solution of the Kosterlitz-Thouless recur- 
sion relations Q 

df(N)/d\nN = -y 2 (N)f(N)/2 
dy(N)/dlnN = [I - f(N)]y(N) , (4) 

where /(iV) = (h 2 nn2D /'2'mkBT)p s (N) / p is a dimension- 
less function proportional to the superfiuid fraction and 
y(N) is proportional to the vortex fugacity. The start- 
ing values fo and yo of the recursion relations ((4]) are 
determined from a best fit to the results correspond- 
ing to different values of N and system temperature 
in close proximity of the transition (see Fig. Q|. From 
these initial values one determines the critical temper- 
ature Tbkt in the thermodynamic limit. Here Tbkt is 
written in units of T* = 27r?i 2 7i2_D/(2mfcB), which pro- 
vides the natural temperature scale for quantum degen- 
eracy in 2D. In Fig. 0] we also show the prediction of 
the universal jump of the superfiuid fraction at the tran- 
sition p s /p = 2mfcsTBKT / '(Trh 2 ri2D) fUj], in nice agree- 
ment with the temperature dependence of our extrapo- 
lated curve. The results for Tbkt as a function of u^dO^d 
are reported in Fig. [5J where they are compared with 
the prediction obtained from Eq. |2]). The agreement is 
surprisingly good up to very large values of the gas pa- 
rameter. At the highest density nir)a& D = 0-1, w hich 
is close to the freezing point ^do^d — 0-32 of a gas of 
hard disks [26] , the dilute gas result ^ is only about 7% 
larger than the PIMC value. 



In conclusion, we have carried out a numerical study of 
the superfiuid transition temperature of interacting Bose 
gases in 3D and 2D and established the limits of validity 
of the asymptotic expansions (JTJ) and @, and the uni- 
versal description in terms of the scattering length. We 
have also explained and resolved the discrepancy between 
previous studies. 
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